Excitation spectrum of a 2D long-range Bose-liquid with a supersymmetry 
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We have studied excitation spectrum of the specfic 2D model of strongly interacting Bose particles 
via mapping of the many-body Schrodinger equation in imaginary time to the classical stochastic 
dynamics. In a broad range of coupling strength a a roton-Iike spectrum is found, with the roton 
gap being extremely small in natural units. A single quantum phase transition between strongly 
correlated supefluid and quantum Berezinsky crystal is found. 



Usually the system of Bose particles at zero tempera- 
ture exists in one of two possible ground-states: super- 
fluid (SF) or crystalline (CR). More exotic option is a 
"supersolid" ground state suggested long ago which 
attracted a lot of attention recently 0; this is a state 
which is expected to possess both superfluid and crys- 
talline order simultaneously. Another direction of the 
search for unusual quantum ground states is related with 
a search for a " Bose- metal" , that is, a bosonic analog of 
a Fermi-liquid, see for example 0, 0] ■ Such a state would 
possess neither superfluid no crystalline order. Sugges- 
tion for the search of such a strange quantum state was 
made 20 years ago in Ref. Q, in relation with classical 
thermodynamics of 3D vortex liquid in high-temperature 
superconductors. This idea was further developed in 
Ref. where two different models of stronglyinteracting 
Bose-liquid were considered (note that Refs. [1,0] refer to 
continuous 2D Bosc-liquids without any lattice, whereas 
Refs. [1,13 consider lattice models). The arguments were 
given in Ref. Q in favor of existence of a new unusual 
ground-state which is still liquid, but is not superfluid. 
One of these models refers to 2D bosons interacting with 
a 2D dynamic U{1) gauge field, with an effective cou- 
pling constant ~ 1. The second model (KKLZ, Ref. 0) 
is purely static, it has a remarkable feature that its ex- 
act ground-state wavefunction is represented in a simple 
Jastrow form. 

It was shown later in Ref. Q that KKLZ model obeys 
nonrclativistic supersymmetry which allows to obtain a 
number of interesting results analytically. The KKLZ 
model contains a coupling constant a such that small 
values a < 3 definitely lead to a gapful superfluid state, 
whereas at very large a > 35 a kind of a "Berezinsky 
crystal" with power-law decay of positional correlations 
is stabilized, according to Ref. Q . An issue was raised in 
Ref. 0] about possible existence of a third, intermediate, 
ground-state of the "normal liquid" type, which could 
exist in some part of the broad range 3 < a < 35. Su- 
persymmetry of KKLZ model makes it also possible to 
compute time-dependent quantum correlation functions 
via classical Langevin dynamics (the relation between su- 
persymmetry and Langevin dynamics was discussed, in 
particular, in Ref. [I!]). Similar approach was proposed 
by C.Henley Uj for the lattice quantum systems and 



used efficiently in Ref. [12[ to explore excitation spectrum 
of quantum dimer models on square and triangular lat- 
tices at the Rokhsar-Kivelson point [l3| . More recently, 
the same lines of ideas were developed in Ref. [l^ for 
quantum spin models. 

In the present Letter we report results of extensive 
numerical studies of dynamic density-density correlation 
function in the KKLZ model through a broad range of 
coupling strength 2 < a < 40. The presence of roton- 
hke branch of the excitation spectrum is demonstrated, 
with the ratio of the roton gap A to the plasma fre- 
quency Wo strongly decreasing with increase of a. Right 
before the crystallization transition at a ~ etc ~ 37, 
this ratio becomes less than 3 • 10^"^; still we could not 
identify any finite interval oi a < ac where roton gap 
would be exactly zero without a Berezinsky crystal be- 
ing formed. An effective roton mass m* defined via the 
spectrum w(p) = A -|- (p — po)^ /2m* near the roton min- 
imum, is found to be weakly-dependent upon a. The 
spectral weight S{p,uj) is well-approximated by the sin- 
gle quasiparticlc peak at w < 2A. whereas at higher ener- 
gies quasiparticle spectrum is undefined due to strongly 
decaying nature of excitations. Our results support an 
existence of superfluid ground-state all the way up to the 
crystallization transition, but the transition temperature 
Tc(a) scales with A (a) and becomes extremely low at a 
close to ac- 

We study the KKLZ model of 2D interacting Bose- 
particles characterized by the exact ground-state wave- 
function of Jastrow form 



^'o(ri, ...,yn) = const 



n 

j>k 



(1) 



Here n is the particle density and a is a parame- 
ter. Many-body probabihty density Po(ri, tat) = 
I'I'olri, rAr)p can be considered as a Gibbs measure 
for a classical 2D liquid with potential energy 
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and temperature T. Quantum Hamiltonian of the KKLZ 
model is defined as 

, a i dV 
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dr 
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where we put T = h and 2m = 1. Langevin dynamics 
leading to the Gibbs distribution Po(ri, r^r) is defined 
as 

dr,,^ _ dV{r,} 
dt ~ 9r,-^ ^' 

where = 2T5jkSf_,^S{t - t'). Our goal is 

to compute dynamic density-density correlation func- 
tion S{k^t) = ^{n\^{t)n_\^{Q)) , (V is the system's vol- 
ume) in the ground state (GS) of the Hamiltonian 
In terms of spectral expansion it is given by 
5(k,t) = ^. |(GS'|nk|k,i)|^e-*'^'^-'* where i denotes ah 
quantum numbers except the momentum k. The equiv- 
alence [1, [13, [m, of quantum and classical dynamics 
for the theories like the one defined by Eq.Q allows us 
to use classical simulation of the Langevin dynamics de- 
fined in Eq.Q to compute S'(k, t) in the imaginary-time 
domain: 5'(k, — ir) = 5(k, r), where 

5(k,r)= j \{dv,{Q)dv,{T) X (5) 

3 

X P(r, (0) , 0, V, (r) .r)^Y. e'"'" (6) 

i,3 

where P(rj(0), 0, rj(T), t) is the 2-time iV-particle joint 
distribution function for the stochastic diffusion process 
defined by Eq.Q. For the derivation of Eq.® see the 
Supplementary material. 

We begin with the application of our computa- 
tional method to the simpler case of the Calogero- 
Sutherland model (GSM) [H] defined on a ID cir- 
cle of the length L. The GSM ground-state wave- 
function is ■00 = Y\i^j^^'^'^{''^Xij / L), where A > 1/2 
and arbitrary otherwise. The corresponding classical 
potential energy is Vcsm = -AX]i<j Insin^Trxi^/L. 
We simulate GSM model with N = 200 particles via 
Langevin dynamics to compute its dynamic structure 
factor S{k, r) and compare with exact results avail- 
able nil- According to Ref. GSM spectral den- 
sity S{k,Lo) = ^J(GS'|Afc|fc, i) 1^ — w) is nonzero 
in a finite region lo G [E^{k), Ej^ik)] only, where 

= Vs {k - 1^) , E+ = Vs (k + , and i;^ = 
ttA^ and k^ = 2TTn. In Fig [T] we plot results of 
numerical simulation for A = 2 together with the- 
oretical low bound curve. In our computation, the 
lower bound of the spectrum was determined as the ex- 
trapolation £'min(fc) = limt^oo din S{k,t)/dt; another 
spectral characteristic is its simple average ujfm{k) ~ 
ujS{uj, k)du! = d\nS{k,t)/dt\t^o- The agreement be- 
tween data for -Emm(fc) and theoretical spectral bound- 
ary E-{K) is remarkable. It proves the capability of 
our method to capture gapless excitations with large 
wavevectors k ^ ko, which are invisible in the "first mo- 
ment" approximation w/m- Note that for small k ko 
data for Emin{k) and ujfm{k) coinside, as it should be for 
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FIG. 1: Computed lower bound of the spectrum i5min(fc) for 
A = 2 CSM is shown by circles, the results for simple spectral 
average iUfm{k) are shown as triangles, exact theoretical lower 
bound (k) is represented by the full line, dashed line shows 
static structure factor S{k). 

the spectral density nearly saturated by single-particle 
excitations. 

Now we turn to our major subject: search for the low- 
energy roton modes in the KKLZ model defined by the 
Hamiltonian ([3]). An example of the excitation spec- 
trum in the strong coupling region, a = 20, is shown 
in Fig. [21 here and below N = 256. We plot here the 
data for E^in{k) for the wavevectors k in the vicinity of 
ko — l-K^/n. where static structure factor S'(fc, t = 0) has 
a peak. The inset to Fig [2] shows E'minlfc) in a broad 
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FIG. 2: Quasiparticle energy near the roton minimum for 
a = 20 obtained by the fit of 5(fc, r); in the inset lower bound 
for -E(fc) in the whole fc-range is shown. 

range of k determined via best fit of S(k, t) to the sin- 
gle exponent Ae~^'-'^. In the main panel of Fig. [2] we 
show E^in{k) in the narrow region around ko, obtained 
via more accurate fifing procedure described in Suppl.2. 

A roton minimum in E'minlfc) is clearly visible at fc = 
ko; below we denote the roton gap as A = i?niin(fco)- For 
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a = 20 the magnitude of the roton gap A is found to 
be very small, about 1% in comparison with the plasma 
frequency ujq = ^Tra^, which sets a natural energy scale 
in the problem. In particular, ujq is the frequency of 
the uniform density oscillations in the KKLZ model, see 
Ref. Q for details. Thus, our first qualitative observa- 
tion is that in the strong-coupling region the excitation 
spectrum shows a very deep roton minimum. As follows 
from the general arguments [17|, a well-defined excita- 
tion spectrum may not exist in the k region where quasi- 
particle decay is allowed by conservation laws. For the 
roton- like spectrum with deep minimum, the " no-decay" 
condition is fulfilled at energies E < 2A only: at higher 
excitation energy, the decay into two rotons is allowed 
with a high rate. A well-defined roton excitations may 
exist in the momentum range p- < k < p+ around the 
minimal point k^. According to Ref. [l7| . the excitation 
energy E{k) is expected to approach the end-points p± 
nonanalitically, with a zero slope: 



E{k^p±) = 2A 



,-6/|fc-p±l 



(7) 



results arc consistent with linear vanishing of the gap at 
a ^ 37 ~ 38, slightly above the point of the crystalliza- 
tion transition aoid = 35 found in Ref. @ for classical 
2D Coulomb gas. However, the values of A in this range 
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FIG. 4: Ratio of the roton gap to plasma frequency in the 
large-Q range. 



where p± are called spectrum terminating points, and a 
and b are some positive constants. The equation ([7|) re- 
sults [TtI from an exact summation of the most singular 
diagrams for the momenta k p±. Our data presented 
in Fig. [5] (main panel) are in good qualitative agreement 
with this prediction; the spectrum end-points arc situ- 
ated at (p„,p+) ss (0.8,1.2) • ko- Unfortunately, high- 
presision computation of i?min(fc) close to the end-points 
was found to the very difficult due to increasing data 
scattering. 

Similar analysis of the relaxation data for different 
vaues of the coupling constant a yields the dependence 
of the gap magnitude A on a presented in Fig. [3] in log- 
arithmic scale. Increase of a leads to very sharp (nearly 
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FIG. 3: Roton gap A as function of a, logarithmic scale 

exponential in the range 10 < a < 35) decrease of the 
gap magnitude A. The same data for the region of large 
a > 20 arc presented in Fig. 0] in linear scale. These 



contain large relative errors which makes it difficult to 
determine unambigously where A(a) vanishes. To ap- 
proach the problem of location of the quantum critical 
point from another perspective, below we compare long- 
time asymptotics of the dynamic structure factor 5(fcoj t) 
in the liquid and crystalline phases. The crystalline 




T, arbitrary units 



FIG. 5: Inverse logarithmic derivative of S{ko,T) is plotted 
for a — 35 (blue dots), 37 (red boxes) and 40 (black bars) as 
measured directly in the simulation. 

phase of the KKLZ model is very specific. This is densly 
packed triangular lattice, but, instead of usual transverse 
phonons with w ~ it supports phonons with parabolic 
dispersion, oj{q) ^ q^. This comes from the fact that 
shear modulus of this lattice vanishes itself in the long- 
wavelength limit, ^{q) oc g^, see Ref. 0; here wavevector 
q = p — Gi, where is one of principal inverse lattice 
vectors. The presense of soft shear modes leads to a spe- 
cific long tail in the time decay of the angle-averaged 
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structure factor S{ko,T) = / ^5(fco cos (p, fco sin r). 
which can be measured by Langcvin dynamics: 



d\nS{kn,T) 1 

7i " o ■l^'^) 

dr 2t 



r » (8) 



^0 



where /(t) > decays exponentially with r and to* is 
the effective mass (to be discussed later). Now we define 
a function y(r) = — (2(iln5(fco, T)/dr)~^ and note that 
according to Eq.(|5]) it should never cross the line y = t. 
On the other hand, in the liquid phase with a nonzero 
gap A, the function j/(r) approaches 1/2A at r — oo, so 
its crossing with the straight line y = t occurs definitely. 
In Fig [5] we present simulation results for the function 
y(r) at a = 35, 37 and 40. According to the criterion 
formulated above, the critical value ac is also found in 
the range 37 < ac < 38. The data summarised in Fig. |4] 
and Fig. [S] support the conclusion that liquid state with 
a small roton gap A transforms into a crystalline state 
via the single phase transition where A vanishes. 

Coming back to the discussion of the the gapful liquid 
phase at a < ac, we note that low- lying excitation with 

« fco are chatacterized, apart from the gap value A, by 
the value of the effective mass to* = (^d'^E{k)/dk^\ko) 
Measurement of the S{k, r) decay in the vicinity of 
ko allows to determine m* in a broad range of a, as 
shown in FiglS] The results shown in Fig. U and Fig. [5] 
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FIG. 6: Quasiparticle mass m* weakly depends on a even at 
the transition point. 

yield the parameters of the low-lying excitation spectrum 
e{k) = A -I- '^^2m* ■ allowing to determine the tempera- 
ture of superfluid-to-normal transition Tc(a). Within the 
Landau-type mean-field theory Tc is defined as the tem- 
perature where superfluid density Us = n — n„ vanishes. 
Neglecting quasiparticles interaction, we find equation for 
the critical temperature Tc'. 



1 
2m 



dfB{e,Tc)^2 d^k 



de 



(2^)2' 



(9) 



valid in the range 10 < a < ac. Note that corrections 
to Tc due to vortex depairing (Berezinsky-Kosterlitz- 
Thouless mechanism) are very weak, due to smallness 
of the roton gap A in comparison with the plasma fre- 
quency OJQ. 

In conclusions, we have computed excitation spectrum 
of 2D Bose-liquid with long-range interaction in a strong- 
coupling regime. Broad range of coupling strengths a < 
ac is found there gapful supcrfiuid state is stable at T = 
in spite of a very small value of the roton gap A. Out 
data suggest a single quantum phase transition from such 
a strongly correlated superfluid into a quantum crystal 
phase at ac « 37 38. At smaller a, superfluid state 
is stable up to the critical temperature Tc ~ 0.3A(a), 
which is orders of magnitude lower than a naive estimate 
To ~ h^n/m would give. 

We arc grateful to L. B. loffe, D. A. Ivanov, L. N. 
Shchur and M. A. Skvortsov for useful discussions and ad- 
vises. This research was supported by the RFBR grant # 
10-02-00554 and by the RAS Program "Quantum physics 
of condensed matter" . 



where /B(e,T) is Bose distribution function. Evaluation 
of the integral (O leads to the result T'c(a) w 0.3A(a) 
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SUPPLEMENTARY ONLINE MATERIAL 

1. Mapping from quantum mechanics to classical 
stochastic evolution 

Standard Fokker-Planck equation corresponding to the 
Langevin dynamics, Eq.(j4]) is 



P(ri 



(10) 

Equilibrium solution of Eq. ([TO| is given hy Pq = e . 
One can check that after the change of variables P = 
^'o(ri, ..rjv)^'(ri, ..Tat, t) the equation (fTO|) assumes the 
form of imaginary-time Shroedinger equation ^' = H'^ , 
where Hamiltonian H is constructed from the potential 
V as shown in Eq.Q. For the following we denote a 
position in coordinate space {yi, ..v^) = <y5 and will not 
use the specific form of H . The correspondence of clas- 
sical and quantum correlation functions that we prove 
below is valid for any symmetric H = which ground 
state "^(^{(f) is known exactly. The symmetry condition 
H = leads to H* = H, which enables us to choose 
real wavefunctions, so P{ip, t) is always real. 

Quantum states form a full system of orthogonal func- 
tions: 



X 

Consider quantum correlation function: 



(11) 
(12) 



Cgit) = i^olAe-'^^'Bl^o) (13) 

Inserting into R.H.S. of Eq. p^ the decomposition of the 
unity operator (|12p we obtain 



/ 



C,{t)=Y,{^o\A\X)e~^''{X\B\^,) ^14) 

A 

A 



where A and B are diagonal operators (i.e. func- 
tions of coordinates tp only). The derivation of the 
quantum-clssical mapping begins with replacing variables 
P{(p,t) = '^o{(p')'^{(p' ,t). The operator governing the 
classical stochastic evolution is W = —'^o{lp)H-^^^. 
It's easy to see that P\{f) = '^^{if')'^ \{ip') are the 
eigenfunctions for this operator, yet this system of eigen- 
functions is neither normalized nor orthogonal since the 
operator is a non-Hermetean one. Combining the 
identity H = and the definition of W , we obtain: 
PqW'^ = WPq, which is the detailed balance condition. 

Rewriting Eg. p^ formally in classical notations, we 
find 

C,{t)= [ dp,d^'Y,Px{v)A{v)Pxy)B{v')e-'^' (15) 



Now we need to evaluate classical correlation function. 
We have the equation for probability density P{(p)'. 



P = WP, 



(16) 



and system is in the equilibrium state P{t) = Pq. The 
two-time correlation function C(t) (as given by R.H.S. 
of Eq.®) is defined via stochastic process 2-time proba- 
bility P{ip,0,ip\T): 

C{t) = / d^V^(^,0,^',T)A(^)i?(^') (17) 



where P((^, 0, v?', r) = P{'p)p[p^^i according to the defi- 
nition of a conditional probability p'^_^^i that the system 
will be in configuration ip' at time r, given that it was in 
configuration (p at time r = 0. Substituting this expres- 
sion for P{(p, 0, p', r) into Eq. pT]) we find: 



C(t) = / d^d^'Po{^)A{p,)pl^.B{^') (18) 



To evaluate p[p^^i — {e^'^)^pi5{ip' — (p), we need to know 
the decomposition of S- function into eigenmodes. it is 
convinicnt to use "quantum" basis (since classical oper- 
ator W is non-Hermetean): 



/\ A 



j:pxi^)Pxi^') 



(19) 



Now we can contract this (5-function with e^"^ (remember 
that eigenvalues are — A) 



EPA(^)PA(^')e-^^ 



^'o(^) 

For classical correlation function we obtain 



(20) 



C(r) = / dipd^'Po{p)A{ip) 



Y.Px{^)Px{p')e 



-B{p>') ^21) 



Po{^) 

^Y.P>.{^)A{^)P^{^')B{p,')e-^^ 



Comparing Eqs. (|2T|) and (fT5|) we find the relation wanted: 
Cq{t)=C{it) (22) 

2. Details of data analysis. 

For rotonic spectrum with gap A quasiparticle contin- 
uum begins at w > 2A. It can be seen by considering 
2 rotons with minimal energy (fci_2 = fco) and arbitrary 
angle between fci and ^2. Total energy is 2 A, and total 
momentum can be set arbitrary in the region K < 2fco. 
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Rotons arc the only detected excitations below the con- 
tinuum: 

Siio, k) = A5{io - E{k)) + Scon{io, k), (23) 
where Sconi^jJ < 2 A) = 0, (24) 

E(k) « A + at k ^ ko (25) 

In the region k G we assume the main contribu- 

tion to come from a quasiparticle, i.e. in Ea.(|23p 

A> y" Scon{^,k)(hj, (26) 

so that the exact shape of S'con(w) does not matter. For 
data fitting we use rectangular spectral density Scon{^) = 
_B, 2A < Lo < uJmax, SO for each value of k there are four 
fiting parameters: A, E(k), B, LOmax, apart from the value 
of A = iJ(fco) that is the same for all k. We minimize the 
mean square deviation fititi) ^ Ssim{ti)Y to find 
E{k) plotted on Fig[5J We also check the condition (^5)) 
and find that it is violated in the close vicinity of ter- 
minating points, thus the statistical error of determining 
E{k) grows there. 

To collect data presented in Fig |31 we do not need to 
use the fc- regions near the terminating points fc « p_ , p+ 
, so we can use inequality (j26p and estimate E{k) just as 
dhiS / dt\t=t„ where to is sufficiently long to lead to ad- 
ditional exponential damping of the continuum modes. 
Note that inaccuracy in determination of ln5 (and of 



its derivative) grows exponentially with to, since 1/iS ~ 
^E{k)ta ^ Therefore the finite simulation time determines 
how long is the optimal interval to we can use. The 
derivative d\nS/dt can be accessed with the use of Monte 
Carlo estimator (subtracting the S values for consequent 
configurations), or by drawing a line through the se- 
quence of points ln5(i,;), t, € [to - ST/2, to + 6T/2] . 
These approaches yield similar results, but the latter is 
more insightful when one tries to assign errorbars aE{k) 
to the results for E{k). These errors contain standard 

N-point slope measurement error a^E = \J~^^^^ and 
the systematical overestimating of E{k) due to the con- 
tinuum modes. The second source if errors is related with 
the fact that \nS{ti) is not exactly linear function of time. 
Assuming that the derivative dlnS/dt changes by SE in 
the interval ST, we can estimate possible systematic er- 
rors as acE - ,5^e-(2A-i5)*T/2/(i _ g-(2A-B)5T)^ The 

denominator of this expression diverges while k approach 
terminating points, which reminds us of the range of ap- 
plicability of the method we used. Surprisingly, the data 
analysis using N-point treatment of IniS(t) and neglect- 
ing systematic shift acE, can be performed in the whole 
range of k € (p_,p_|_). This method catches non-analitic 
behaviour of E{k) near the spectrum terminating points, 
as well vanishing of the roton gap A while a approaches 
ac- In both these cases, errors bars unE grow consider- 
ably, indicating the approach to a transition. 



